Systematic experimental exploration of bifurcations with non-invasive control 
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We present a general method for systematically investigating the dynamics and bifurcations of a physical 
nonlinear experiment. In particular, we show how the odd-number limitation inherent in a popular non-invasive 
control scheme, (Pyragas) time-delayed feedback control, can be overcome for experiments with periodic 
forcing. To demonstrate the use of our non-invasive control, we trace out experimentally the resonance surface 
of a nonlinear oscillator near the onset of instability, around two saddle-node bifurcations (folds) and a cusp 
bifurcation. 
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Feedback control is not only of interest as a tool for system 
manipulation in the classical control engineering sense, but 
it can also be used for model verification and discovery, if 
one can ensure that it is non-invasive. Pyragas' time-delayed 
feedback (TDF) is a popular feedback control scheme that 
is automatically non-invasive [1]. It feeds a signal u{t) = 
k'^\y{t) —y{t — T)] back into the experimental dynamical sys- 
tem, where y is some system output (possibly processed), k is 
a vector or matrix of control gains and T is an a-priori chosen 
time delay. If the time delay equals the period of a periodic 
orbit ^(f) (f e [0,r]) of the uncontrolled dynamical system 
and the gains are such that the control is stabilizing, then the 
controlled system will also have the periodic orbit y, because 
the control input u vanishes for all time (that is, the control 
becomes non-invasive). However, the control has changed the 
stability of y, making it visible in the experiment. 

While sometimes TDF (or its extended version [2]) is used 
for engineering purposes (suppression of period doublings 
leading to chaos [3]), non-invasiveness is not essential in these 
applications. TDF draws interest mostly in the scientific com- 
munity because it enables experimenters to explore dynamical 
phenomena such as equilibria and periodic orbits of the origi- 
nal uncontrolled system regardless of their dynamical stability 
[4, 5]. 

Systematic studies that try to explore parameter space of the 
uncontrolled system and that use TDF to stabilize equilibria 
and periodic orbits non-invasively encounter a major difficulty: 
it is not known under which conditions one can find control 
gains k that successfully stabilize a periodic orbit [6]. This 
is in contrast to classical feedback control where one feeds 
u{t) = k^\y{t) —y^,{t)] with a pre-determined reference signal 
yH.(f) back into the system. For classical feedback control it 
is known that, if corresponds to an equilibrium or periodic 
orbit of the uncontrolled system then, under some genericity 
assumptions (controllability and observability), one can always 
locally stabilize the equilibrium or periodic orbit even if it has 
arbitrarily many unstable eigenvalues [7]. Experimental and 
theoretical studies have explored the limits of applicability of 



TDF and restrictions on the gains due to instabilities caused by 
the TDF feedback term F \y{t) - y{t - T)] [5, 8]. In particular, 
there are topological restrictions (the odd-number limitation 
[9]) which guarantee that TDF cannot possibly work in some 
cases. Unfortunately, one common scenario where it would 
be natural to use TDF is ruled out by the odd-number limita- 
tion: an equilibrium of an autonomous system or a periodic 
orbit of a periodically forced system in the vicinity of a system 
parameter setting where it makes a fold (called saddle-node 
bifurcation, see the black markers in Fig. 2 for a typical bi- 
furcation diagram). We note that even the unstable controller 
proposed in [10] fails to stabilize uniformly near the fold. 

In this paper, we demonstrate a simple alternative ap- 
proach to explore bifurcation scenarios, including the unstable 
branches. The demonstration experiment is the forced nonlin- 
ear oscillator shown in Fig. 1, an electro-magnetic energy har- 
vester [11, 12], mounted on a force-controlled electro-dynamic 
shaker. 
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FIG. 1. (a) A schematic of the physical parts of the nonlinear energy 
harvester, (b) A schematic of the experimental set-up. The elements 
within the shaded box are implemented within a real-time control 
system. See [11] for implementation details. 
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FIG. 2. Experimental data showing the evolution of the controlled 
system as the bifurcation diagram of the uncontrolled system (shown 
in Fig. 1(a)). A family of periodic orbits is tracked through two saddle- 
node bifurcations (folds). Artificially large steps are taken along the 
solution curve for illustration purposes. Forcing frequency: 22 Hz. 



Our approach exploits the fact that the goal of the experiment 
is a parameter study in a system parameter p (in the set-up of 
Fig. 1, p{t) — acos{cot), so we have two parameters), rather 
than finding a single equilibrium or periodic orbit at a specified 
parameter value. It is applicable whenever the feedback control 
is achieved by (effectively) varying the same system parameter 
p that one wants to use for the bifurcation diagram. We explain 
the basic principle first for equilibria, as illustrated by the inset 
(a) in Fig. 2 (the underlying data comes from the experiment, 
which is a harmonically forced oscillator, and requires a small 
modification of the basic principle, discussed afterwards). 

Suppose we have an experiment with a single system param- 
eter p (which we can vary) and a single output y, and we know 
that the feedback control law u{t) — +k[y^: —y] (control gain 
A; 7^ 0) is stabilizing the output at parameter value for all 
{p*,yt) along a branch of equilibria. Assume also that we have 
found already two previous points along the branch of equilib- 
ria, {pn-i,yn-i) and {pn,yn) (labeled by circular markers in 
Fig. 2(a)). We set 

(PreiJief) = {Pn,yn) + h[{pn,yn) " {Pn-Uy„-l)] (1) 

(h = 1 in Fig. 2(a)), and run the experiment with control law 

M(f) =/'ref + ^bref-3'(f)]- (2) 

The experiment will settle asymptotically to a constant output 
value jasy, which we call yn+i, and a constant input value 
^asy — Pvef + ^bref — Jasy], which wc Call Pn+i, such that the 
new point along the branch is (/7„+i ,y„+i ). Only the combined 
parameter Fref — /^ref/^ + 3'ref is varied during the experiment 
at discrete times. After each change of Yy^f and some transient 
time, the asymptotic values of the inputs and outputs give the 
parameters p„+i and the output values yn+i. In Fig. 2 two 
phases of the transients are visible between the black markers: 
the horizontal coordinate p can be set, so there is an initial 
sharp increase of p, followed by a gradual stabilization due to 
the control toward {pn+i,yr:+i)- 
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FIG. 3. Time profile of forcing and response amplitudes and error 
(non-harmonic part of control input ii). Snapshots are time profiles 
corresponding to inset (a) of Fig. 2. Note that the time gaps between 
the time profiles are only gaps in the time series recordings due to the 
saving of data (typically of the order of milli-seconds); the experiment 
ran continuously. Forcing frequency: 22 Hz, sampling frequency: 
5 kHz. 



The approach illustrated in Fig. 2 should be compared to the 
general approach of solving the nonlinear fixed point problem 



Jref =}'asy(/'ref,3'ref), 



(3) 



in combination with a pseudo-arclength condition using a New- 
ton iteration, introduced in [13]. The brackets in (3) indicate 
that the asymptotic output jasy depends locally uniquely on 
parameter and reference signal if the feedback is stabilizing. 

The method promises a substantial speed-up compared to 
[12, 13] by removing one equation per free bifurcation parame- 
ter from the fixed point problem, instead of adding an equation 
as is typically the case. The projection onto the solution surface 
occurs along the line {{p,y) '■ p — fref + ^bref — y]} determined 
by the control gains k. Whenever the remaining equations of 
(3) can be solved by a fixed-point iteration (or there are no 
equations remaining), this removes the need for a full Newton 
iteration. 

When the basic idea above is adapted to the scenario of a 
harmonically forced oscillator (where the bifurcation parameter 
is the forcing amplitude) the combined parameter Y = p + ky,:ef 
becomes time-dependent (periodic with the forcing period 
2n/co). One can still apply feedback law (2) using (/^refj^ref) 
defined by (1), and the input u{t) will converge to a function 
Masy(f) asymptotically. However, Masy will not be harmonic 
because the asymptotic output 3'asy(0 is not harmonic for a 
nonlinear oscillator. In this case, we can remove the higher 
harmonics using an iteration: define P/^x as the ^th (complex) 
Fourier coefficient of a periodic signal x, such that x{t) — 
£^^oRe[P<;Xe''^™]. Then we iterate the reference signals y,ef = 
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FIG. 4. Experimental results from the energy harvester shown in Fig. 1 . A sequence of constant forcing frequency runs were performed at a 
spacing of 0.2 Hz. Panel (a) shows the complete resonance surface of the oscillator. Panel (b) shows the corresponding two-parameter bifurcation 
diagram (a top-down view of panel (a)) with the cusp point evident at approximately 19.2 Hz. Panel (c) is a front view of the resonance surface 
with the measured error superimposed onto the surface. The error is defined as the root-mean-square (RMS) difference between u{t) and p{t) as 
a percentage of the forcing amplitude; it measures how invasive the method is. In all panels, the data points are shown as black dots and the 
calculated saddle-node bifiu'cation (fold) curve is marked in black. Points within the dark gray region of panels (a) and (b) are unstable solutions. 



Pief/k + ymf in a Picard iteration: 

Y,,f^o{t)^Y„it)+h[Y„it)-Y„^i{t)] (4) 

Fref,j+1 (f ) = Re[Pl (FrefJ -yasyjOe'™] +y.syj{t), (5) 

where y^sy.j{t) is the asymptotic output signal for reference 
signal Frefj CO- That is, the non-harmonic part of the next 
reference signal Fref is set to the non-harmonic part of the 
asymptotic signal of the output jasy. 7(f)' ^^d the harmonic part 
is kept constant for iterations j > I- 

In the experiment shown in Fig. 1, the non-harmonic parts 

of 

^asy C3,n be removed (to experimental accuracy) in a single 
iteration for all parameter values considered (see Fig. 3 and 
Fig. 4(c)). Thus, the overall protocol to find points along the 
solution branch is as follows. 

1. Set yref,o(f) = Y„it)+h[Y„it)-Y„_iit)]. 

2. Run the experiment with control law u{t) — PD(3'(r) — 
Fref,o(0) until transients have settled, then record the 
asymptotic output 3'asy,o(f )■ 

3. Set Fref.l (f) = Re[Pl (}^-ef,0 -3'asy,o)e'™(f)] +3'asy,0(f 

4. Run the experiment with control law u{t) ^ PD(3'(f ) ~ 
iref,i(0) until transients have settled, then record the 
asymptotic input Masy.i (f ) and output yasy.i (f )■ 

The next point on the branch is then (OiJn+i (0) = 

("asy,i(0i3'asy,i(0)' where pn+\{t) is harmonic up to exper- 
imental accuracy and yasy,i(0 i^ asymptotic output of the 
experiment in step 4. Note that the Fourier decomposition of 
u and y does not need to be done in real-time and instead can 



be done as a post-processing step to choose the new reference 
signal Fref.i and to check convergence. 

Figure 3 shows the time series recordings corresponding 
to the data points shown in Fig. 2(a); they demonstrate in 
detail the convergence of the method as the system passes 
through a saddle-node bifurcation (fold). The first harmonic of 
the input Piu (Fig. 3(a)) gradually drifts during non-periodic 
transients but rapidly settles. The error e{t) — u{t) — Pi Me"*" 
(Fig. 3(b)) corresponds to the non-harmonic, invasive, part of 
the control; its decay during each step is evident. The vertical 
bars (marked fo, ti and t2) indicate the stages of the iteration: 
step 2 occurs from fo to t\, and step 4 occurs from ty to t2- The 
input u and output y at t2 are then accepted as points on the 
branch when their corresponding Fourier coefficients become 
stationary for 5 consecutive forcing periods. Figure 3(c) shows 
the amplitude of the first harmonic Pyx of the displacement to 
further demonstrate convergence. 

The main advantage of the method presented here over meth- 
ods based on Newton iterations, apart from ease of implemen- 
tation, is the speed-up of a factor of k, 15 compared to [12, 13] 
(a conservative estimate; only individual solution curves could 
be traced out in [12, 13]). This feature is particularly important 
if one wants to explore systems that gradually degrade under 
laboratory conditions. 

As illustrated in Fig. 4(a), this speed-up enables tracking 
of entire surfaces and the associated bifurcations. The exper- 
imental data points (marked by black dots in panels (a) and 
(b)) are obtained by consecutive runs for fixed frequencies 
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0.2 Hz apart. The total experimental time to generate these 
results was 61 minutes. Panel (a) shows the three-dimensional 
projection in the space spanned by the two parameters forcing 
frequency and amplitude and the response amplitude (note 
that the response is non-harmonic and so this is indeed only 
a projection). Its main feature is the curve of saddle-node bi- 
furcations passing through a cusp bifurcation (black) [14, 15]. 
Curves of constant forcing amplitude (gray), reminiscent of 
the resonance curves for an idealized Duffing oscillator, give 
additional geometric information. All the data points within 
the dark gray shaded area of Fig. 4(a) are unstable periodic 
orbits of the uncontrolled system, and as such would typically 
not be seen experimentally. 

Figure 4(b) shows a top-down view of panel (a), a two- 
parameter bifurcation diagram, again indicating all measured 
points on the unstable part of the surface in a darker shade 
of gray. The saddle-node bifurcation (fold) curve bounds the 
instability region with a cusp point at approximately 19.2 Hz. 

Figure 4(c) shows a front view of panel (a) with the error 
at each data point rendered onto the surface. Here the error 
is defined as the root-mean-square (RMS) difference between 
the actual input signal u{t) (as fed into the shaker) and the 
reported harmonic forcing signal p{t) (the harmonic part of 
m) as a percentage of the forcing amplitude. In essence, this 
is a measure of the invasiveness of the control; if this method 
was truly non-invasive, then this error would be zero. In the 
experimental set-up here the error is low, with a mean error of 
< 0.5%. The predominant source of error is noise amplification 
through the use of a derivative controller. This error is kept to 
a minimum through the use of an appropriate digital filter [11]. 
As seen in Fig. 4(c), there is no apparent correlation between 
geometric features of the solution surface (e.g., the fold points) 
and the magnitude of the error at that point. 

The presented approach is a general method to explore dy- 
namical systems in parameter studies near saddle-node bifur- 
cations (folds). It is particularly useful for the exploration of 
families of equilibria because no iterations similar to (5) are 
necessary. As we demonstrated, it is also applicable to peri- 
odically forced systems (the generalization to a non-harmonic 
forcing is straightforward). The main limitation of the method 
is that control fails at the linear level whenever the system does 
not depend on the bifurcation parameter to first order (e.g., near 
transcritical bifurcations). As the presented approach works 
particularly well around saddle-nodes, its main application 
areas are likely complementary to those of Pyragas control. 
Examples currently under investigation include the identifica- 
tion of growth rates in chemostats [16], or tracking localized 
spots in ferrofluids [17] 
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